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ABSTRACT 



A continuous identification of parameters is performed 
on a simulated fast breeder nuclear reactor system using 
hybrid computation and applying techniques of statistical 
regression analysis and exponentially -mapped-past functions. 
Output states which are not directly measurable are estimated 
by use of a Kalman filter. 

The method developed in this study is applied to a nu- 
merical example which demonstrates that unknown parameters 
can be identified within 3% of their actual value, with sig- 
nal noise ratios as lov>f as 10:1 in the measured states. The 
example also demonstrates that convergence occurs in a rea- 
sonably short time. 

" An executive software routine was written in FORTRAN IV 
for the XDS-9300 digital computer and was applied both for 
the control of the process and for establishing the opera- 
tional performance and accuracy of the Comcor Ci-5000 analog 
computer that simulated the nuclear reactor system. 
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I. INTRODUCTION 



Fast breeder nuclear reactors such as the recently 
proposed design described in [1] do not have inherent load 
following ability and therefore require an automatic con- 
troller to handle load changes. It is important in design- 
ing for stable operation of the automatic control system to 
know the value of the parameters of the plant within a rea- 
sonable degree of accuracy. These parameters change due to 
such factors as temperature, fuel depletion, fuel loading 
and power level, [2], making it desirable to have continuous 
identification of parameters. 

The continuous identification of parameters of dynamic 
systems has received intensive attention during the last 
few years [3], [4]. References [5] and [6] among others, 

have shown successful results by applying gradient techniques 
and quasilinearization, respectively. These methods are cum- 
bersome for nonlinear systems and in particular for dealing 
with a nonlinear space -dependent reactor system, [6], where 
the requirement of continuous identification cannot be met 
since the in-core sensors which are necessary for certain 
state measurements have not yet been developed^. 

The implementation for continuous parameter identifica- 
tion is difficult and complex. For example, many gradient 

1 

The procedure available consists of the exposure of 
a matallic foil by immersion and takes enough time (between 
one and two hours) to eliminate any possibility for con- 
tinuous identification of parameters. 
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techniques, even when they represent in some cases very ele- 
gant .ways of approaching the hill-climbing search procedure 
are limited by practical considerations in their application 
To illustrate this, reference is made to a particular gradi- 
ent (steepest descent) technique which uses an error crite- 
rion and sensitivity equations, and which, among other 
requirements, has the need for a forcing function. A block 
diagram of the algorithm is shown in Figure 1. In this 
method it is required to derive, implement and solve the 
sensitivity equations in order to obtain the value of the 
gradient for the specified cost function. In order to iden- 
tify n parameters the number of analog components, using 
standard methods, must be multiplied by at least n * , since 
for a linear case these equations take the same form of the 
model, [7]. In addition, for nonlinear systems the sensi- 
tivity equations have more complicated forms which do not 
permit the use of techniques such as multiplexing, [7]. To 
apply these methods the system must be initialized with some 
previous information, usually an educated guess of the pa- 
rameter values, and then the process must be iteratively 
reinitialized until some criterion is satisfied. Since the 
parameter correction is a function of the gradient and 
since this gradient approaches zero as the parameter value 
approaches the actual value (solution) , the procedure will 
converge in a reasonable time only if a means is provided 



More recent techniques, [8], allow for simplification 
when the model can be expressed in canonical form. 
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Figure 1. Identification of parameters by gradient technique using sensitivity 
equations . 



to speed up the convergence in the vicinity of the minimum 
of the cost function. 

Other techniques, [9] and [10], which are used in non- 
linear systems suffer from more or less the same problems. 
These difficulties indicate the need for an identification 
system that will operation without supervision, require 
neither educated guesses nor initialization, be very reli- 
able, and use simple implementation that avoids the trouble 
some iterative initializing process. 

The procedure developed in this thesis makes use of 
the method of statistical regression analysis, applied pre- 
viously with success to linear systems, [4], and confined 
for a long time to the analysis of statistical populations, 
[11] and [12] . The method of steepest descent is applied 
to the classical regression analysis in order to speed up 
the convergence. To avoid the need for reinitialization, 
the properties of exponentially-mapped=past statistical 
variables is used, [13] . 

Although the procedure is applied here to a simplified 
model of a fast breeder nuclear reactor system, the method 
is general enough to be applicable to many other automatic 
control systems that are represented by more complex models 

For the particular nuclear reactor system used as the 
example application, the variable representing the delayed 
neutron concentration is a nonmeasurable state and in order 
to have its value for the identification system a Kalman 
filter, [14] and [15], was used. 
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The unknown parameters were identified within 3% of 
their actual values. The example demonstrates that con- 
vergence occurs and that the time of convergence compares 
favorably \vith current methods, [16]. The simulation re- 
sults shown in this thesis were obtained using zero as the 
initial values for the parameters; however, tests have 
sho\m that the convergence time can be greatly improved if 
an educated guess is accepted. 

An executive software routine was written in FORTRAN 
IV for the XDS-9300 digital computer and was applied both 
for the control of the process and for establishing the 
operational performance and accuracy of the Comcor-5000 
analog computer that simulated the nuclear reactor system 
and was used to implement the analog part of the identifi- 
cation procedure. 
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II. REACTOR DYNAMICS 



The point model of the nuclear reactor kinetic equa- 

2 3 5 

tions for U fuel can be found in several references, 
[16]-[19]. 

These equations are 



and 



dn 6K 



B 






n. I 



i=l 



\.C. 
1 1 



(1) 



dCi B^ 

“dTt 1 = 1 , •••, 6 . 



( 2 ) 



The left-hand term in equation (1) represents the rate 
of change of the neutron density; the first term on the 
right side of equation (1) represents the rate of production 
of prompt neutrons, the second term the rate of absortion 
plus leakage and the last one, the production rate of de- 
layed neutrons. The variable 6K, the reactivity input, can 
be expressed as 



6K = p 






j-1 



.T. . 
1 1 



Applying this substitution, equations (1) and (2) become 

1 6 



dt 1 ^e 



Z oc.T, I ^ - 4 n + Z X.C, 



j=l 



'i 3 



i=l 



1 1 



C3) 



dC. 

1 

“at 



^i 



n - X.C. 
1 1 



(4) 
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where T is defined by 



dT. / 

j ^ ^ j ( K7 • ^ 




j=l 



( 5 ) 



Since and are fractional quantities it follows from 
their definitions that 



and 





= 3 . 



The initial conditions for equations (3) thru (5) are 



n(0) 
C. (0) 



n 



C. 

lO 



^i % 

Z* X. 
1 



i = 1, 



,6 



T.(0) = 0 



1 = J- 






For the purpose of this study a single- time-constant 
feedback reactor with a single group of delayed neutron 
precursors will be used; the following equations define the 
system 



dn 

3t ■ 





n 3 
' Z^ 



n + XC 



( 6 ) 



dt Z* 



XC 



(7) 



e 





oT. 



( 8 ) 
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The initial conditions are 



p;(o) = 0 

n(0) = n 

e n 

= Co r*ir 

T(0) = 0. 

By changing variables as indicated in [16] , (see Ap- 
pendix A) a more suitable representation is 



N = ^ [(Kg - aT - 1) N + D] 


(9) 


D = X (N - D) 


(10) 


T = Y (N - 1 - bT) , 


(11) 



where 



N 






6x10^ 




Y = 



y 

e 




and the initial conditions are 



N(0) = 1 T(0) = 0 

D(0] = 1 Kg(O) = 0. 

Experience with present day design of nuclear reactors 
indicates that these systems are stable under normal opera- 
ting conditions. The large breeder reactors are still in 
the design stage and a computer simulation analysis has 
shown that these plants will be stable, [1]. Nonetheless, 
a stability analysis applying the method of Liapunov has 
been performed for this model and is included in Appendix B. 
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The model defined by equations (9) thru (11) was im- 
plemented, but acceptable accuracy could not be obtained. 

The large value of 3/Jl* required a large scaling factor and 
also a large gain for the state variable N. Because of 
these two reasons the voltage that represented the state 
variable N become too small and too close to the uncertainty 
level (fourth digit) of the analog computer. 

From simulation it is observed that N is approximately 
equal to zero. This and the fact that is very large 

implies that equation ( 12 ) is a satisfactory approximation 
for equation (9). This further simplification \\rill, of 
course, introduce some inaccuracies since it is assumed, in 
fact, that N is zero; however, it has been shown, [16], that 
the largest error introduced with this assumption is less 
than 21 . 

With this simplification the model is 



N = - D/(K 



e 



aT - 1) 



( 12 ) 



D = (N - D) 



(13) 



T = T(n - 1) - byT. 



(14) 



In state variable form where 



X 



D 



Xi = X (X 3 - xj 

^2 = Y (X 3 - 1) - byx 2 



(15) 



(16) 



and the variable X 3 is defined as 



X3 = - Xj /(K^ - ax2 - 



ax 2 



1 ). 



(17) 
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For a numerical example, the known constants are as- 
sumed to be k = 0.4 and y = 1.0 and the parameters to be 
identified are "a” and ”b". Using these values the equa- 
tions are 



Xj = 0.4 C ^3 - Xj) 



Xj = Xj - 1 - bX2 



( 18 ) 

(19) 



and the variable Xj is defined in equation (17) . The analog 
computer equations using the scaled variables 



and 



Xj = 40xi 
X 2 = 40x^ 

X 3 = 2x 3 



are then 

Xj = 8 Xj - 0.4x^ 

X 2 = 20Xj - 40 - bx^ 
where the variable x is defined by 



( 20 ) 

( 21 ) 



X 



3 



2Xi 

40K - ax 2 - 40 

e 2 



( 22 ) 



Figure 2 shows the implementation of the above equa- 
tions and the settings of the potentionmeters for this 
scaling are shown in Table I. 
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Figure 2. Analog implementation of the nuclear reactor simulation. 



TABLE I 



POTENTIOMETER SETTINGS 



Pot 


Value 


Pot 


Value 


pool 


0.200 


P012 


0.200 


P002 


0.800 


P015 


0.600 


P003 


0.500 


P017 


0.600 


P004 


0.500 


P021 


0.600 


POOS 


0.200 


P051 


0.600 


P006 


0.400 


P055 


0.600 


POlO 


0.200 


P057 


0.600 


poll 


0.600 
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III. REGRESSION ANALYSIS 



Several attempts were made in this study to achieve 
the desired goals for the identification system as stated 
in Section I. Arrising in all of these attempts ivas the 
problem stated by R. C. K. Lee, [14], and illustrated in 
Figure 3, which is related to the closed-loop formed by 
the estimator and parameter identification system. This 
problem is difficult because the overall system is nonlin- 
ear even for linear plants. Because convergence of this 
loop was not proved, and because estimation and identifica- 
tion of parameters for nonlinear cases depends strongly on 
the initial value^, it was decided to approach the problem 
in the following alternative way. 

The system is nonlinear as defined by equations (9) 
thru (11). However, the simplification performed on the 
model that leads to equations (18) and (19) yields a second- 
order linear model with the variable Xj representing an in- 
put. Since it is possible to measure N and T, which are 
represented in the simplified model by the variable Xj and 
the state variable x^ , respectively, then the non-measurable 
state D, represented by x^ , can be easily estimated. In 
order to make the simulation realistic, the variable Xj is 
treated as a non- deterministic forcing function. This func- 
tion is generated by applying equation (22) and the resulting 

2 

It is possible for this kind of system to reach a 
steady-state value different from the actual value (this is 
known as a state of nonlinear stability) . 
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value is corrupted with some Gaussian noise o£ mean zero 
obtained from the low-frequency Gaussian noise generator 
in the Comcor-5000 analog computer. 

This simplification has the advantage that it opens 
the estimator-identification loop, Figure 4, because the 
only estimated state (x^) is neither a function of the un- 
known parameter "a”, nor a function of the unknown parameter 
"b”. 

The technique used in this study, known as statistical 
regression analysis, consists basically of finding the best 
fit to a set of data using a leas t- squares criterion. It 
has been shown, [4] and [10], how the steepest descent 
method can be used on the analog computer for parameter op- 
timization and model building, and therefore this method 
was applied in the derivation of the regression equations 
(27) and (30) . 

In order to perform the regression analysis it is re- 
quired to define S, a cost function. The cost function, 
for simplicity, is usually chosen to be the integral of 
the square of the error, which is defined by a relationship 
that contains the parameter to be identified. The second 
step is to find the rate of change of S with respect to a 
new independent variable^, x, and to minimize the cost func- 
tion using the method of steepest descent. The details of 
these methods are developed in the following sections. 



The variable x represents computer time and it changes 
more rapidly than real time. 
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A. IDENTIFICATION PROCEDURE FOR THE UNKNOWN PARAMETER B 

.From equation (14) an error function Cj can be defined 
as 



Ej - T - Y(N-l-bT) 



and a cost function Sj as 



Si = 



= J' ^t = J' [T-T(N-l-bT) ] 2 dt 
0 0 



(23) 



(24) 



is an estimate of the parameter b. 

The partial derivative of Si with respect to time is 

3Sl _ 9S, 

8x 8b 8 t 



Using the method of steepest descent applied directly 
to the objective function it is required, [4], to make 



3t ab 

/\ 

This ensures that as t increases the estimate b changes 
in such a way that Si decreases. When the estimate b makes 
3b/9T = 0 the minimum of the cost function has been reached. 

Taking the partial de.rivative of equation (24) with 
respect to b gives 



3Sj 

— t : — 

3b 

and, therefore 




[T- Y(N-l-bT) ]y T dt 
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9t 



t t t t " 

^Ttdt-rC NTdt + Y I Tdt+b f T dt . (25) 



Ttdt- Y NTdt + Y 




Letting 



X 



D 



X 



2 



T 



(26) 



X, = N 



and substituting (26) into equation (25) gives the final ex- 
pression 



Since x^ and x^ are both measurable and continuously 
available, the identification of the unknown parameter "b” 
could also be done continuously and the input to the iden- 
tification system does not require any digital processing. 

Due to a lack of analog multipliers, however, all multipli- 
cations were performed in the digital computer. Figure 5 
shows the implementation diagram for equation (27) . In 
Section III-D a modification to the procedure is introduced 
in order to avoid the need for reinitialization; this modi- 
fication uses the properties of exponentially-mapped-pas t 
averages. 

B. IDENTIFICATION PROCEDURE FOR THE UNKNOWN PARAMETER A 

The parameter "a” can be identified using equation (12) . 
Since D is a non-measurable state and it is needed in the 



3t 




(27) 
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Figure 5. Identification system for parameter "b" using regression analysis and 
steepest descent only. 



regression equation, a Kalman filter, as described in Sec- 
tion III-C, was used to obtain the estimate of the variable 

D. 

Taking equation (12) 

N = - D/CKg-aT-1) 

and defining an error criterion as 

£2 ^ D + N(K^-aT-l) , (28) 

leads to the cost function 

t 

^2 = / ^2 ^ 29 ) 

0 

a is an estimate of the parameter "a". The steepest descent 
algorithm can be applied to obtain the differential equation 
that defines the regression analysis for the unknown parame- 
ter "a". 

Using (26) gives 






and 






9a 9 a 



t 

I 



+ X, (K^-a X, - 1) ] ^ dt 






Therefore , 
9 t 



t 



= 2 






L 



X2 X 3 dt + J X 2 X3 dt 



28 



t t 




( 30 ) 



0 0 



Figure 6 shows the analog implementation of equation (30) . 

As in the identification of parameter "b”, (Section 
III-A) , this implementation was augmented by using the 
properties of exponentially-mapped-pas t variables (see 
Section III-D). This was done in order to avoid the trouble 
some iterative reinitialization. 

C. KALMAN FILTER 

As stated in Section II, the linear equations (15) and 
(16) are used in the Kalman filter for the estimation of the 
state variable D. These equations are repeated here for 
convenience . 



In the application of the Kalman filter the variable Xj is 
treated as a forcing function. 

Figure 7 shows the block diagram of the estimator for 
the system defined by the difference equations 



where k represents a point in the discrete time domain, _x(k) 
is the vector of the state variables at time k, w(k) is the 
system input which in this case is the sum of two components 



Xi = A(X3 - Xi) 



(15) 



X2 = ■ 1 ■ bX2) . 



(16) 



x(k=l) = $ X (k) + r w (k) 
^(k) = H X (k) + V (k) 
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Figure 7. Block diagram of the Kalman filter. 



the variable Xj (a deterministic input) and Gaussian noise. 
z,(k) is the vector of measured states at time k, and v(k) 
is the measurement noise. The matrices , T and H are, 
respectively, the state transition matrix, the distribution 
matrix and the measurement matrix. 

Recurrence relations for the optimal filter, [14], 
[15], are 

x(k|k) = £(k|k-l)+G(k) [z(k)-H x(k|k-l)] 
x(k+l |k) = J x(klk) , 

/s 

where x(k|k) is the estimated state vector at time k using 

/N 

observations up to and including time k, ^(k|k-l) is the 
predicted state vector at time k, using observations up to 
time k-1 and G(k) is the gain matrix. 

The gain matrix G(k) is given by 



^(k) = ^(k|k-l)H^[H P(k|k-1)H^ + R] 



- 1 



P(k|k-1) is the covariance matrix of prediction error defined 
by 

^(k|k-l) ^ E{[x(k) - x(k|k-l)][x(k) - x(klk-l)]'^}. 

The auxiliary equations that determine the propagation 
of P are 

£(k]k) = (I^ - £(k) H) _P(k|k-l) 

P(k+l|k) = f P(k|k) 4 ”^ + Q 

where JP (k I k) , called the covariance matrix of estimation error, 
is defined by 
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P(k|k) = E{[x(k) - ^(k|k) ] [x(k) - x(k|k)]^}. 



R, the covariance matrix of the measurement noise, and Q, 
the covariance matrix of the system input noise, are de- 
fined as 



R = E[v(k) v"^ (k) ] 

0, = Jf E[w(k) w’^(k) ]^ . 



The matrix form of the difference equations, that cor- 
respond to the differential equations (15) and (16) , used 
for the filter, is 



x(k+l) = 



• 




m 4 






.992 0’ 




.008 0l 




X3 




x(k) + 








0 .990 




o 

• 

o 


1 


X -1' 


■ m 






1 


3 ^ 



H = 



0 

1 



for a sampling interval of .02 sec. and a value of .5 for 
the parameter "b”. 

Since "b" is required for the $ matrix then the iden- 

/s 

tification of "a" is done after the identification of b has 
been completed. 

The signal- to-noise ratio is specified as the ratio of 
two exponentially-mapped-pas t averages, the square of the 
signal and the square of the noise signal. This average 
function is explained in the next section. 
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The computer listing of the estimator algorithm con- 
taining the complete set of equations that characterize the 
estimator is within the subroutine ESTIM. 

D, THE EXPONENT IALLY-MA.PPED- PAST AVERAGE FUNCTION 

The open- loop approach used for the implementation of 
the Kalman estimator produced system equations which do not 
depend explicitly on the unknown parameter "a”. This imple- 
mentation of the Kalman estimator allows use of the regres- 
sion analysis technique. However, this procedure allows 
the variable Xj to reach the identifier without being fil- 
tered. Another problem still present is the iterative re- 
initialization. 

To solve these problems an exponentially-mapped-past 
(EMP) average function, [13], v\ras found very useful. The 
EMP function also was found to be very suitable for acting 
both as a filter and eliminating the overloading of the in- 
tegrators used in the identification. This exponential 
weighting function decreases with early values of time and 
its implementation allows for relatively fast control of 
its time constant. It gives the advantages of easy com- 
putation and simple implementation. 

The EMP average function is defined as 



where u(t) is a unit step function. This convolution in- 
tegral of f(t] and an impulse response of a low-pass RC 



00 




- oo 
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filter, with leakage time constant RC = 1/a is shown im- 
plemented in Figure 8. 

Figures 9 and 10 illustrate the changes made in order 
to include the EMP function into the simulation. The 
changes are the addition of input potentiometers with set- 
tings of a. 
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Figure 8. Exponentially-mapped-past average 
circuit . 
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Figure 9. Identification system for parameter 
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Figure 10. Identification system for parameter " 



IV. EXECUTIVE ROUTINE 



One of the most important tasks prior to obtaining 
data from the analog simulation is the verification of the 
operational performance of the analog computer. This is 
especially critical where the analog computer used in the 
simulation is large and has many other users, [20]. This 
task is very difficult and no attempt was made here to cover 
the whole system; however, it was found worthwhile to check 
out all components used in this simulation. 

The executive routine MAIN can be divided into two 
subsections: the first one performs the checkout and the 

second one controls the actual identification. 

The test performed by the checkout routines are: 

1. TESTO. - This subprogram tests the library sub- 
routines SETPOT, SCAN and DAC . To accomplish this TESTO 
sets potentiometers and compares the desired setting with 
the actual setting, allowing for the inherent uncertainty 
present in the fourth digit. In the same way the digital- 
to-analog converters (DAC's) are tested for proper setting 
and also for any interaction among them. Any detected 
error will cause a message to be printed and a PAUSE to be 
executed . 

2. TESTl. - This subporgram tests the initial con- 
dition setting on the plant and the analog-digital conver- 
ters and digital-analog converters. It takes 1000 samples 
of each set of converted values and prints the mean errors 
and the error variances of both tests. 
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3. TEST2. - This subroutine tests the integrators 
used in the simulated nuclear reactor system. It prints 
out the percentage error of the integrator integration 
rates . 

4. TEST3. - This subprogram tests the initial condi- 
tions on the simulated plant, which allows verification of 
the divider M003. It also runs the plant to compare the 
final value with a precomputed numerical solution. This 
subroutine prints the percentage error for both tests. 

5. TEST4. - This subroutine performs a test of both 
identification systems using a precomputed numerical solu- 
tion for the steady-state and prints the percentage errors 
of both systems compared with the parameter value used in 
the numerical solution. 

The identification of the unknown parameters is per- 
formed by subroutines RUNl and RUN2. 

1. RUNl. - Activates the identification system for 
the parameter ''b'' , which is performed first because the 
Kalman filter requires the value of the unknown parameter 
"b". The internal subroutine SAMPl, activated by interrupt 
53, performs the sampling of X 2 and X 3 , computes the deriva 
tive of x^ by a difference approximation, sets the values 
of the reactivity input and the required products in the 
simulated nuclear plant and identification system respec- 
tively. RUNl prints on the teletypewriter the value of "b” 

2. RUN2. - Activates the identification system for 
the parameter "a". HANDY, an internal subroutine activated 
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by interrupt 53, samples X 2 and X3, calls ESTIM (Kalman 
filter) to obtain an estimate of xi, performs the required 
multiplications and sets these values in the analog computer 
for the identification system. 

The program was written primarily in FORTRAN IV for 
the XDS-9300; some assembler language statements were used 
in order to be able to use the library subroutine DELAY. 

A sample of the output taken under conditions of 
failute of the library hybrid executive routine is shown 
as a part of the Computer Output. 



41 



V. SIMULATION RESULTS 



Several tests were performed for each identification 
system. The first test of each system is shown in Fig- 
ure 11 and Figure 12. These figures show the response of 
the identification system to a constant input whose value 
was numerically precomputed and corresponds to a steady- 
state value of the plant, i.e., the digital multipliers 
shown in Figures 9 and 10 were furnished with these steady- 
state values. Since these runs were performed without state 
variable trajectories they do not reflect the previously 
specified desired objective of an identification system. 
Nonetheless, for this example the identification system 
converges to the actual parameter value if supplied with a 

s 

noiseless discrete sample representing the steady-state 
value of the plant variables. The time scales are un- 
labeled since the convergence time can be altered arbi- 
trarily in any noiseless situation by adjusting the system 
gain, without changing the final result. Figure 15 shows 
the type of plant transient condition used, a truncated 
ramp. For the identification of the unknown parameter "b” 
a low gain was used and for the parameter "a” a high gain 
was used. These choices of lovj and high gain were made 
arbitrarily because for noiseless runs the procedure was 
not sensitive to gain and either high or low gains produced 
the same satisfactory results. Here again, for the reason 
stated previously, there is no need for labeling the time 
scales . 
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Figures 16 thru 23 show the behavior of the system 
under simulated real conditions. 

The system response in the identification of the parame- 
ter "b” with two different levels of measurement noise is 
shovm in Figure 16, where the noise was added to the vari- 
able that represents the plant temperature, . 

Figures 17 and 18 show the system response in the iden- 
tification of parameter "a.” with two different levels of 
measurement noise. Unbiased noise was added to the vari- 
able that represents the plant temperature, , and to the 
variable that represents the neutron density, x^. 

Figures 19 thru 21 v^^ere obtained with biased measure- 
ment noise, and Figures 22 and 23 show typical inputs to 
the identification system. 

Talbe II summarizes the numerical results for various 
runs with different values of the parameters "a” and "b”. 
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TABLE II 



NUMERICAL RESULTS 



PARAMETER IDENTIFICATION 



a 


b 




a 


b 


0.4 


0.3 




0.40 


0.29 


0.4 


0.4 




0.40 


0.39 


0.4 


0.5 




0.40 


0.50 


0.4 


0.6 




0.41 


0.61 


0.4 


0.7 




0.40 


0.71 


0.4 


0.8 




0.39 


0.81 


0.2 


0.5 




0.21 


0.50 


0.3 


0.5 


0 


0.29 


0.50 


0.5 


0.5 




0.50 


0.49 


0.6 


0.5 




0.59 


0.50 


0.7 


0.5 




0.69 


0.50 



The signal- to-noise ratios used in all these runs was 20:1. 
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Figure 11. Identification of parameter "b" under steady - 
state conditions. 
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Figure 12. Identification of parameter "a" under steady- 
state conditions. 
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Figure 13. Identification of parameter "b” under transient 
conditions . 







Figure 14. 



Identification of parameter "a” under transient 
conditions . 
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Figure 15. Reactivity function used in the simulation. 
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b) Signal to Noise Ratio 10:1 

Figure 16. Identification of parameter ”b" under noisy 
conditions . 
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Figure 17. Identification of parameter "a” with measurement 
noise on . 
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a) Signal to Noise Ratio 15:1 




b) Signal to Noise Ratio 10:1 



Figure 18. Identification of parameter "a" ^^^ith measurement 
noise on X 2 . 
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Figure 19. Identification of parameter "b" under a biased 
noise condition. 
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Figure 20. Extreme biased condition for parameter 
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Figure 21. Extreme biased condition for parameter 
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Figure 22. Typical input signals for the identifier of 
parameter ''b”. 
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Figure 23. Typical input signals for the identifier of 
parameter "a”. 
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VI. CONCLUSIONS 



From the results of the numerical example used, the 
following conclusions are reached regarding the desired 
goals for the identification procedure. 

A typical response of the identification system where 
the initial value of the parameters are set to zero is shown 
in Figures 16 thru 21. The results indicate that the pro- 
cedure does not require a guessed value to be effective. 
Other runs were made with various initial parameter values 
and if the guessed value was close to the actual value, the 
convergence time was greatly reduced as might be expected. 

In any case, whether or not a non zero initial value was 
used, the system did not require iterative reinitialization. 

The simplicity of the identification procedure imple- 
mentation, together with the small number of required analog 
components, provides for a very reliable system, since ana- 
log integrators and analog multipliers are very reliable 
components . 

The average convergence time for a signal- to-noise 
ratio of 10:1 was found to be approximately 16 seconds, and 
about 10 seconds for noiseless measurements. If the signal- 
to-noise ratio is small the convergence time is large. 

For the example considered in this study the procedure 
identified the unknown parameters within 3^ of their actual 
values with signal- to-noise ratios as low as 10:1 in the 
measured states. 



56 



The identification system for parameter ''b'' was very 
insensitive to the noise, except for biased noise (with a 
mean of 10 % of the stead-state value of the state variable) , 
as shown in Figure 20. The extreme biased condition was 
taken as the maximum noise for which, even when the system 
does not actually converge, the identifier output oscillates 
within 4% of the actual parameter value. 

The identification of parameter "a" also showed very 
good performance but it was more affected by measurement 
noise on X 3 than on X 2 . This was expected since X 3 is not 
filtered by the Kalman filter and X 3 is an important factor 
in every term of equation (30) . This problem can be elim- 
inated by using filtering in the measuring circuits for X3. 
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APPENDIX A: DYNAMIC EQUATIONS SIMPLIFICATION 



As seen in the case of an heterogeneous reactor made 
up of n independent media, the state variables are related 
by the following non linear differential equations: 




n 

I* 



^ n 



Vx . c . 



( 1 . 1 ) 



dC. 3 . 

1 « 1 ^ n p 

dt Z* ^ ‘ ^i^i 



( 1 . 2 ) 



1) -t -,(T. - 



T^). (1.3) 

The initial conditions for this set of equations are 



n(0) 
C. (0) 
T.(0) 



n 



10 

0 . 



l*X 



n 



For the purpose of this study the complexity of the 
model is reduced to a single time constant approximation, 
i.e., i = 1 and j = 1, so that 



^ = (Pe - “T) . XC (2.1) 

ar'lln-^C (2.2) 



- oT 



^ dt _ / n 
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(2.3) 



with the following initial conditions 



P^CO) = 0 
n(0) = n^ 



‘ vS: 

T(0) = 0 . 



A further manipulation will produce a more convenient 
form, i . e . , 



dn _ _3 

3T a* 



_3 

I* 



or 




Jl ^ 

n^ dt I* 3 

In the same way 




(3.1) 



dc 



^ n - AC 




n - 






= A 





and 



1 






e 



dt 

3t 





- oT 



^ _ U / n 

dt c in 
' o 





(3.2) 



(3.3) 
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Now if 



N 



A _n 



n 




C 




then 

N = ^ [(K^-aT-1) N + d) 

D = X (N - D) 

T = Y (N-l-bT) , 

and 



N(0) = 1 T(0) 

D(0) = 1 K^(0) 



0 

0 . 



A P 

K = -T 
e 3 



, A a 

b = — , 

y ’ 



C4.1) 

(4.2) 

(4.3) 



60 



APPENDIX B: REACTOR STABILITY 



As it was shown by Liapunov, "Problems general de la 
stabilite du mouvement," the general problem of stability 
of motion relative to certain quantities can be reduced in 
many cases to a solution of problems on stability of the 
trivial solution of the system. In theorem 1, "Direct 
Method," he stated: "If the system of differential equa- 

tions that defines the system of differential equations 
that defines the system is such that it is possible to find 
a positive definite function V, the derivative of which, 
with respect to time, calculated by virtue of the system 
equations satisfies the inequality dV/dt < 0, then the sys- 
tem is stable." 

The problem is to investigate whether or not the 
stationary state is stable with respect to initial distur- 
bances . 

The operation of the reactor (mathematical model) is 
defined by equations (3.1), (3.2) and (3.3), assuming that 
the external control is disconnected and neglecting the 
poisoning effects (since they represent a very slow phenom- 
enon) the reactivity can be considered as some function of 
temperature, i.e.: 



6K = - f (t) . 

First of all it should be noted that n and C are always 
positive (and so are N and D) , since they cannot be negative 
for physical reasons. Now if the initial conditions are 
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n = n(t ) > 0 

C = C(t ) >0 
o - 

T = T(t ) , 

then the solutions n(t) and C(t) which correspond to these 
initial conditions satisfies the inequalities 

n(t) > 0 
C(t) > 0 

so the stability will be investigated in the domain D(n > 0, 
C > 0) , inherent stability in the presence of a single group 
of delayed neutrons 



dn 

dt 


-f(T) n - n + AC 


(5.1) 


dC 

dt “ 


n - AC 


(5.2) 


dt 

dt " 


T(n-l) - xbT 


(5.3) 



which obviously admit the special solution for the stationary 
state 

T = Tj some final temperature (6.1) 

n = n^ = 1 + bT (6.2) 

C = C^. (6.3) 

Now if 3>0, £*>0, X>0, b>0, 1+b =0 are verified then 
this solution for the system is asymptotically stable; for 
physical reasons the first four inequalities of the statement 



62 



are always fulfilled and the fifth is a consequence of the 
fourth one. 

Considering the Liapunov function 

T n C 



V = 



f(t)dt + 



T(n-n ) 
n 



dn + 



T(C-C„) 



dc (7) 



in which the constants T, , n and C are given by formulas 
(6.1), (6.2) and (6.3), the function V becomes zero for 
T = T,, n = n and C = C and under the conditions given in 
Theorem (1) by M. V. Popov in [5] it is positive in the 
domain n > 0 , C > 0 . 

From equations (5.1), (5.2) and (5.3) the total deriva- 
tive of V with respect to time is 



dV _ 
It ~ 



(n-n ) 

f(t) [r(n-l)-rbT] • r ^ ° 
* r f p- n - XC) 



[-f(T)n-j| n + XC] 



= -rbTf(T) - ^ [ IT n - XC] [C(n-n^) -n(C-C^) ) 

= -rbTf(T) - ^ n - XC (8) 

and again under the conditions of the same theorem, the 
function above is negative semi-definite within the range 
n > 0, C > 0 and become zero for 

js" ri = XC 

where T = T i . 
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Since dV/dt is negative semi-definite, it can be con- 
cluded that the solution is stable in the sense of Liapu- 
nov . 
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COMPUTER OUTPUT 



The computer output shown here corresponds to different 
runs where failures in the XDS-9300 library hybrid executive 
package were detected. 
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ERROR IN P057 
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NEUTR9N 

DELAYED 
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FINAL VALUE PERCENTAGE ERRSRS: 

DELAYED NEUTRONS CONCENTRATION 
TEMPERATURE 

neutrons concentration 



testa 

identification percentage errors: 
parameter b 
parameter a 



.4A188 

. 28^60 

.29163 



10.00000 

2.50000 
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CeMPUTE percentage ERR9R AND PRINT MESSAGE 




November 29, 1972 



U.S. DEPARTMENT OF COMMERCE 
National Technical Information Service 

5285 Port Royal Road 
Springfield, Virginia 22151 



United States Naval Postgraduate School 72-5 

Monterey, California 93940 



Re AD 738 874 Solution of a Nuclear Reactor Parameter Identification 
Problem 

Gentlemen: 

We are unable to provide the public with satisfactory copy of referenced 
document for the reason(s) checked below. Please help us by sending 
replacement copy together with this letter by AIRMAIL to address shown 
below. If you are unable to furnish replacement, please return this 
form by AIRMAIL indicating where the necessary replacement copy is 
obtainable. Your cooperation will be appreciated. 



l^Ti 1. These pages are missing, please furnish: Page 75. 

□ 2. These pages are not reproducible, please send good pages: 

□ 3. Document is not reproducible, please send good copy. 

□ These pages have been mutilated, please send replacements: 

□ 5. Other 
Input Section 

National Technical Information Service 



To: Mr, Thott 

I regret that our Library copy 
is also incomplete in that it lacks 
page 75, too. We have contacted the 
advisor in a search for the original 
but it appears that it has either 
been discarded or was retained by the 
author, who is an officer in the 
Peruvian Navy and who haS now returned 
to Ms ooontr,. ^ 

"Georg^ R, Luckett 
Librarian 



Department of Commerce 
5285 Port Royal Road 
Springfield, Virginia 22151 

S^^erely, ^ > 

& 10 ' 

b6 w. thott 

Chief, Input Section 
703-321-8517 
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SUBROUTINE TO CONTROL THE IDENTIFICATION PROCESS 
FOR THE PARAMETER B 
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2$ CALL SETLINES(A, 1.0) 
TL20=0- 1 

CALL SCAN( AHA013>3) 



CALL HSLD 
B = 100»^»B 
OUTPUT(lOl) 
RETURN 
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SUBR0UTINE SAMPl 
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CALL 0AC(1/ .2000/2/ .4000/ 12/ - .4000) 
CALL RESET(IOOO) 

1=5000 

CALL Cer^PUTE 
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